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Abstract 


Quantitative measures of the uncertainty of Earth System estimates can 
be as important as the estimates themselves. Second moments of estimation 
errors are described by the covariance matrix, whose direct calculation is im- 
practical when the number of degrees of freedom of the system state is large. 
Ensemble and reduced-state approaches to prediction and data assimilation 
replace full estimation error covariance matrices by low-rank approximation- 
s. The appropriateness of such approximations depends on the spectrum of 
the full error covariance matrix, whose calculation is also often impractical. 
Here we examine the situation where the error covariance is a linear transfor- 
mation of a forcing error covariance. We use operator norms and adjoints to 
relate the appropriateness of low-rank representations to the conditioning of 
this transformation. The analysis is used to investigate low-rank representa- 
tions of the steady-state response to random forcing of an idealized discrete- 
time dynamical system. 


Correspondence to: ff 


1 



1 Introduction 


The value of an estimate of the Earth system state is increased when accompanied by a 
measure of its uncertainty. Sophisticated users of weather and climate predictions base 
decisions on both forecast and forecast uncertainty (Changnon et al., 1999). Data assimi- 
lation systems combine information from different sources in a manner that depends upon 
its presumed uncertainty (Cohn, 1997). Ensemble prediction systems generate ensembles 
of initial conditions with statistics that attempt to reflect analysis errors (Barkmeijer et al., 
1998). 

Uncertainty can be modeled as a random variable which is completely described by a 
probability density function (pdf). A partial description of the pdf is given by its mean 
and covariance. With some simplifying assumptions, equations can be obtained for the 
evolution of the mean and covariance. However, for realistic dynamical systems with a 
large number of degrees of freedom, the direct solution of these equations, particularly 
those for the covariance, is impractical due to computational costs and to incomplete 
knowledge of the sources of uncertainty. Therefore, in many situations the covariance 
must be modeled. 

One approach to covariance modeling is to specify analytically parameters such as the 
variances and correlation lengths (Dee and da Silva, 1998; Rabier et al., 1998). Complex 
features such as flow dependence may also be modeled through appropriate parameteri- 
zation (Riishpjgaard, 1998). A second approach is to assume that the error variability is 
described well by just a few dominant structures and hence that the covariance matrix is 
approximated by a low-rank matrix. Low-dimensional covariance representations are of- 
ten directly connected to the dynamics, as in reduced-state Kalman filter data assimilation 
and ensemble prediction (Cane et al., 1996; Cohn and Todling, 1996; Houtekamer and 
Mitchell, 1998; Molteni et al., 1996). 

Comparison of low-dimensional covariance representations with the full covariance is 
only possible in idealized models where the number of degrees of freedom is small and 
the sources of uncertainty are specified (Cohn and Todling, 1996; Kleeman and Moore, 
1997; Whitaker and Sardeshmukh, 1998). Here our approach is to examine the equation 
satisfied by the covariance and to determine which of its properties control the covariance 
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spectrum. We consider the situation where the error covariance is a linear transformation 
of a forcing error covariance and seek the properties of the linear transformation that affect 
the covariance spectrum. 

The paper is organized as follows. Section 2 introduces low-rank representations and 
the linear covariance equation. Section 3 analyzes this equation using first the singu- 
lar value decomposition and then operator norms. Section 4 applies the analysis to the 
discrete-time algebraic Lyapunov equation and illustrates the results with an example us- 
ing the dynamics of a generalized advection equation. Conclusions are given in Section 

5. 

2 Mathematical setting 

2.1 Covariance 

We suppose that our system state e, a real vector of length n, is a mean-zero random 
variable. In typical problems, e might be forecast or analysis error, or climatological 
anomaly. Considerable information about the system state is contained in the n x n 
covariance matrix P defined by 

P = (ee T ) (1) 

where () T denotes transpose and (•) denotes expectation. The total variance of e is giv- 
en by tr P where tr denotes trace. The eigenvectors of P order state-space directions 
according to the variance they explain 

A measure of how well a rank-fc matrix approximates the covariance matrix P is 
the quantity |]P — P^||; || • || is a matrix norm, for instance the Schatten p-norm defined 
by 

/ n \ l /P 

|X||, = (X>?(X)j ■ 1<P<°°. (2) 

where < 7 j(X) is the i-th singular value of X. This family of matrix norms includes several 
common matrix norms; for covariance matrices, ||P|| t = trP and ||P||oo = Ax(P) where 
Aj(P) is the i-th eigenvalue of P ordered in decreasing modulus. 
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The optimal rank-A; approximation P*A) of P in the Schatten p-norm is its projection 
onto its leading k eigenvectors given by 

k 

P [k) = '52*i(P)w i wT (3) 

i=l 

where Wi is the orthonormal eigenvector of P corresponding to the eigenvalue A*(P). In 
this case, the error of the approximation is 

n 1 /P 

E A ?( p ) > w 

i— fc+1 

and the appropriateness of low-rank representations of P depends on its spectrum. How- 
ever, direct calculation of the eigenvalues of P is often not practical. 

2.2 Linear covariance equation 

Our approach is to to identify properties of the equation satisfied by P that control the 
spectrum of P. We consider the situation where the covariance matrix P satisfies an 
equation of the form 

P = CP 0 ; (5) 

£ is a linear operator on matrices and P 0 is a forcing covariance matrix. Many systems 
where the covariance evolution is linear and the forcing is stationary can be written in 
this form. An example of a system with this form comes from the stochastically-forced 
discrete-time linear system 

Ck + i = AfcCfc + b k+ 1 , {bibj) = SijP o ; (6) 

e fc is the system state at time-step k, A k is the dynamics, b k is the random forcing and % 
is the Kronecker delta. The covariance P* defined by 

P fc = {e k e T k ) - (7) 

satisfies the discrete-time covariance evolution equation 

k 

P fc+ i = AfcPfcAfc + P 0 = Y, A*Po(A fe ) T = £P 0 , (8) 

i= o 
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which has the form of (5); we use the notation A fc = AkAk-i • • ■ A^. 

In the special case of stable 1 time-independent dynamics A k = A, the equation for the 
steady-state covariance matrix P is 

£- x P = P — APA t = P 0 . (9) 

We will examine in some detail equation (9) known as the discrete-time algebraic Lya- 
punov equation (DALE) (Kwon et al., 1996; Tippett and Marchesin, 1999a, b; Tippett 
et al., 2000a). We refer to £a as the DALE operator; this nomenclature and notation 
is differs from that of Tippett et al. (2000a). Covariance evolution in a continuous-time 
system can also be expressed in the form of (5) (Byers and Nash, 1987). 

We next examine what properties of (5) cause its solution P to have a good low-rank 
approximation. Since the initial or forcing covariance P 0 is often poorly known, results 
that require knowing details of P 0 are unsatisfactory. Results that depend mostly on £ 
and require limited knowledge of P 0 are desirable. 

3 Analysis of linear matrix equations 
3.1 Inner products 

Linear algebra methods such as the eigenvalue and singular value decompositions can be 
used to analyze (5), treating the covariance matrices as vectors of length n 2 and £ as an 
n 2 x n 2 matrix (Byers and Nash, 1987; Ghavimi and Laub, 1995). Eigenmode analysis 
is appropriate when the linear operator £ is normal , that is when £ has a complete set of 
orthogonal eigenvectors. The normality of £ depends on the definition of orthogonality 
and hence on the choice of inner product on matrices. Any inner product (•, •) on vectors 
can be written in the form 

(^i> e 2 ) = ejMe2 (10) 

where M is a Hermitian positive definite matrix and () t denotes conjugate transpose. The 
adjoint A* of A with respect to this inner product is A* = M -1 AtM. We may assume 
’The discrete-time dynamics A is stable if and only if the eigenvalues of A lie inside the unit circle. 


5 



M = I without loss of generality since a new state variable e = M 1/,2 e can be introduced 
with (e 1; e 2 ) = e\e 2 . 

A natural matrix inner product, though not the most general, is defined by 

(X,Y) =trX t MY. (11) 

The matrix inner product (11) is compatible with the vector inner product (10) in the 
sense that the orthogonality of two rank-1 matrices eie\ and e 2 e 2 is equivalent to the 
orthogonality of ei and e 2 . The DALE operator £a has the pleasant properties with respect 
to this inner product that C* A = £a* and that £a being normal is equivalent to A being 
normal. We mention finally that for M = I the matrix inner product (1 1) is the Euclidean 
inner product on vectors of length n 2 since 

n n 

trX+Y = J^XyY a- (12) 

i — 1 j = 1 

3.2 Eigenvectors and singular vectors 

The singular value decomposition of a linear operator gives its rank, range, null space, 
2-norm and optimal low-rank approximations. Approximations of £ that are low-rank in 
the space of nxn matrices are obtained from its singular value decomposition; when £ is 
normal it is sufficient to use the eigenvalue decomposition. Suppose that U;, and V* 
are respectively the i-th singular value, left singular vector and right singular vector of £. 
That is, the n x n matrices U* and \£ satisfy £\£ = cr i (£)U i and (Uj, Uj) = (V, : , Vj) = <% 
where (-, •) is the matrix inner product. Then the decomposition of P in the left singular 
vectors of £ is 

n 2 

P = £P 0 = ^ °i(£) (Vj, P 0 ) Uj . (13) 

i= 1 

Approximate solutions of (5) are found using low-rank approximations of £. If £ has 
some large singular values, for instance, suppose that cr fc (£) » cr fc+1 (£), then a natural 
approximation P of P is obtained using the optimal (in the space of nxn matrices) rank-A; 
approximation £ (k ^ of £ to obtain 

k 

P = £ {k) Po = Y. a ^ ( V - R o) U, - (14) 

i==l 
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with error 


||P-P|| 2 = £P„= CT i(£) l(Vi, Po) I • (15) 

• i=fc+l 


So if the projection of P 0 onto the first k right singular vectors of £ is not too small, 
P shares with P a large component in the directions spanned by the first k left singular 
vectors of £. However, there is no guarantee that the approximation P is low-rank unless 
there is some information about the rank of the left singular vectors of £. For instance, 
if the singular vectors of £ are rank-1 matrices, then P is at most rank -k, though not, in 
general, the optimal mnk-k approximation P (fc) . In general the singular vectors of £ are 
not rank-1 matrices. 

An example where £ is normal and its eigenvectors are rank-1 matrices is the normal 
DALE operator £ A . The eigenvalues and eigenvectors of £ A are (1 — X i (A)X j (A))~ 1 and 
ZiZj respectively where A* (A) and z { are respectively the eigenvalues and eigenvectors of 
A (Lancaster, 1970). The DALE operator £ A is normal if and only if A is normal, in 
which case the eigenvectors z { of A form a complete and orthonormal basis, and P 0 can 
be decomposed in this basis as 

n n 

p 0 = X! (4 p o*i) Ziz) . (16) 

i=l j = 1 


The representation 




i j = i 


Ai(A)A,-(A) 


Z iZ ] 


(17) 


of the solution of (9) shows that there are two extreme situations where the covariance 
matrix P has a dominant low-rank part; in one case the spectrum of P is determined by 
that of A, in the other by P 0 . 

In the first situation A is nearly unstable with some eigenvalues near the unit circle. 
In this case, P has large components that project onto the leading eigenvectors of A, if 
the projection of P 0 onto the leading eigenvectors of A does not happen to be small. For 
example, for P 0 = I and normal A, the solution of (9) is 


p = £ 


1 

1 — |Aj(A)| 2 






(18) 
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and the eigenvectors of A are eigenvectors of P. In this case, P has a good low-rank 
representation if A has a few, but not all of its eigenvalues near the unit circle. The second 
situation occurs when P 0 has some large eigenvalues and the spectrum of A is fairly flat. 
Then P has large components on the space spanned by the leading eigenvectors of Po- For 
instance, when A = cl,0<c<l,P given by 

p = T—— £ p 0 > (19) 

has a good low-rank representation if P 0 does. These two mechanisms can interfere with 
each other and examples can be constructed where P 0 has large eigenvalues and A is 
nearly unstable but the spectrum of P is flat. 

There is no correspondingly simple analysis of the DALE for nonnormal dynamics. 
However, in the next section we show that features seen in the analysis of the normal 
DALE operator hold both for the nonnormal DALE operator and for general £. Namely 
for the DALE we show the relation between the stability of A and the spectrum of P and 
for general £ the connection between the conditioning of £ and the spectrum of P. 

3.3 Operator norms and positive maps 

Examining the eigenvalue and singular value decompositions of the linear operator £ 
does not use the fact that £ is an operator on matrices and hence does not necessarily 
give information about how £ relates the matrix properties of Po and P. This kind of 
information can be found in the operator norm of £ defined by 

l|£!l ' = “ ax W' (20) 

The operator norm ||£|| p relates the matrix properties of P and P 0 , though, for general 
operators there is not a simple expression for its value. 

A linear operator £ with the property that it maps covariance matrices to covariance 
matrices is positive map and hence has the property that its oo-norm is simple to charac- 
terize and that its conditioning is related to the spectrum of P (Bhatia, 1997; Tippett et ah, 
2000a). Many properties of solutions of the DALE and of the continuous time algebraic 
Lyapunov equation (CALE) are direct consequences of the solutions being given by pos- 
itive maps (Bhatia, 1997). For instance, the fact that positive maps have their maximum 
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response in the p = oo norm to uncorrelated forcing 


l|£|l “ = m pf7^y = ra pf W^) =Xl(a) ’ (21> 

has consequences for solutions of the DALE and CALE (Hewer and Kenney, 1988; 
Gahinet et al., 1990; Kenney and Hewer, 1990). The eigenvectors of the “bound ma- 
trix” B = £1 order the state space directions according to the possible response of the 
system there and can be used in estimates (Tippett and Marchesin, 1999a). 

Similarly, results relating the forcing P 0 that produces maximum variance to the re- 
sponse of the adjoint system to uncorrelated error forcing (Kenney and Hewer, 1995; 
Farrell and Ioannou, 1996; Kleeman and Moore, 1997) also hold for any positive map. 
Since the 1-norm measuring variance response is dual to the oo-norm, ||jC||i = ||£*||oo = 
||£*l||oo. Additionally, the linear operator £ achieves its 1-norm on the rank-1 matrix 
ww T where w is the leading eigenvector of £*l since 

||£n;n; T || 1 = (l, Cww T ) = w T (C*\)w = Ai(£*l) . (22) 

The eigenvectors pf the B r = £*l are the stochastic optimals, ordering state space 
according to how much variance is excited by forcing there as shown by the relation 
tr P = tr£P 0 — (Po, Bt) (Bhatia, 1997). 

The operator norm of £ is also related to low -rank representations of P. A necessary 
condition for P to admit a low-rank representation is that it be ill-conditioned, that is, for 
the condition number «oo(P) to be large where ^(P) = |(P|| 0 o||P~ 1 l|oo The conditioning 
of P depends on that of £ and of Po as shown by the relation 

Koo(P) < «oo(£)Koo(Po) ; (23) 

more general results are available for £ A (Tippett et al., 2000a). Therefore, for the special 
case P 0 = I, P can be ill-conditioned and have a good low-rank approximation only if £ 
is ill-conditioned. 

Another connection between low-rank representations of P and the operator norm of £ 
is seen in the relation 

Ai(P)/trP lM 

MPoVWp, * IWU|C 11 ' <24) 
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which follows directly from the operator norm definition and bounds the fraction of the 
total variance explained by the first eigenmode of P. A significant fraction of the total 
variance of P can be in its first eigenmode if the same is true of P 0 , or if the quantity 
||£|| 00 ||£~ 1 ||i is large. 

An additional useful property of positive maps is that they preserve ordering 2 . That is, 
if one has upper and lower bounds Pq < Po < Po > then 

£Pq < P < £Pj . ' (25) 

These bounds give bounds for the eigenvalues, diagonal and total variance of P 

Ai(£Po ) < Ai(P) < A*(£P+) , (26) 

diag(£Po ) < diag(P) < diag(£Pjj-) , (27) 

tr£Pg < tr P < tr^CPo - . (28) 

Using the upper and lower bounds with Pq = A n (P 0 )l and Pq = Aj(P 0 )l gives bounds 
that depend on the bound matrix: 

A n (Po)B < P < A 1 (P 0 )B. (29) 

When the bounds in (29) are tight and P has a well-separated set of leading eigenvalues, 
the leading eigenvectors of P and B span approximately the same subspaces (Golub and 
Van Loan, 1996, Theorem 7.2.4). However, if P does not have a well-separated set of 
leading eigenvalues, P - A 1 (P 0 )B may be small without the leading eigenvectors of P and 
B spanning approximately the same subspaces. 

4 Analysis of the DALE 

The preceding analysis can be used to determine when the solution of the DALE (9) 
admits a low-rank representation. The condition number of the DALE operator requires 
calculating the norm of £ A and For the normal DALE operator, ||£ A || p = Ax(£ a ) = 
(1 — |Ai(A)| 2 ) -1 and (Tippett et al., 2000a). 

1 - |A n (A)| 2 < ||£ a 1 Up < 1 + | A x (A) j 2 , (30) 

2 For two symmetric matrices X and Y, X < Y means that Y - X is positive semidefinite. 
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showing that £ A is ill-conditioned when A has some but not all of its eigenvalues near the 
unit circle. 

In the case of nonnormal dynamics A, ||£ A 1 || p can be estimated in terms of the singular 
values of A, 

max |1 - of (A) | < W^Wp < 1 + A) (31) 

and ||£ A || P can be bounded by 

2r(A) + r 2 (A) “ ,|£a|,p ~ r(A)2 (32) 

where the radius of stability r( A) is the distance from A to the closest unstable matrix 
(Mori, 1990; Tippett et al., 2000a) 3 . These estimates show that the DALE operator £ A is 
ill-conditioned when A is nearly unstable and the A has at least one singular value that is 
not near unity. 

An example of stable nonnormal dynamics that produces a highly ill-conditioned DALE 
is the dynamics coming from the generalized one-dimensional advection equation 

Ht + a/Ar + c'(x)n = 0 , 0 < x < 1 (33) 

with a > 0 and initial and boundary conditions 

= 0) = Hq{x), (i(x = 0, t) == 0 , (34) 

respectively. Nonnormality is due to the undifferentiated term c'{x)^jl and the boundary 
conditions. A similar model with periodic boundary conditions is described in detail in 
(Tippett et al., 2000b). Disturbances move from left to right with speed a. As disturbances 
pass through regions where d{x) < 0 they grow and where c'(x) > 0 they decay. The 
zero boundary condition at the left boundary forces the solution to be identically zero 
after the time o _1 required to cross the domain. As a consequence, the operator A T that 
advances the solution r time units is nilpotent and all its eigenvalues are identically zero; 
there is no modal growth. 

3 For normal stable matrices r(A) .= 1 — |Ai(A)|, the distance from its largest eigenvalue to the unit 
circle. For nonnormal dynamics the radius of stability depends on the pseudospectrum of A (Trefethen, 
1997).Eigenvalues near the unit circle, large singular values and sensitive eigenvalues cause the radius of 
stability to be small. 
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The discrete-time dynamics A r that advances the solution r time units is given by (see 
Appendix) 

i 

{ 0 0 < x < ar 

(35) 

s(x)/j,(x — ar, t) ax < x < 1 , 

where s(x) = exp [(c(x — ar) — c(x))/a]\ note that s(x) = exp— rc'(x) in the limit 
ar <C 1. In this example we take a = l/12,r = l and c(x) given by 

1 2 

c(x) = 1 + — — arctan(16(0.5 — x)) . (36) 

16 7 r 

In most parts of the domain c'(x) « 0 so that there is little growth. Around x — 0.5, 
c'(x) < 0 and there is amplification as shown in the plots of c(x) and s(x) in Fig. 1. 

Nonmodal transient growth is found from the singular values and singular vectors of A T . 
The singular values and left singular vectors of A T are the square roots of the eigenvalues 
and the eigenvectors of A r A^ (see Appendix) 


AA T ^(x, t ) = 


{• 

(^s 2 (r)/^(r,f) 


0 < x < ar 
ar < x < 1 . 


(37) 


The singular values of A are zero and the values taken on by s(x). The maximum growth 
in one time step is cri(A) = 1.37 is sensitive to c'{x). Where s(z) > 1 there is local 
amplification. In the spatially discrete case, the singular vectors of A are columns of the 
identity matrix. The leading singular vectors are associated with the region where s(x) 
has its maximum. 

We add mean-zero Gaussian distributed random noise bk at each time-step: 


Cfc+i — A T ek T bk , (pb T ) — Po 


(38) 


with covariance P 0 = 1 + G where G is a Gaussian correlation model with correlation 
length 0.25 and normalized so that tr G == 1. A sense of the time behavior of the system 
is seen by looking at the mean with respect to x of the forcing bk and of the response y k in 
Fig. 2 which show that the dynamics amplify the forcing and increase the time coherence. 

The maximum possible amplification is found from the diagonal matrices B = £1 and 
Bt = £at\ to give that ||£a||oo = 1209.6 and ||£a||i = 1157.6. The norm of £ A is 
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large because though the eigenvalues of A are far from the unit circle, they are sensitive to 
perturbations. The maximum of B is within ar of the right boundary and the maximum of 
B t occurs within ar of the left boundary, in contrast to the leading singular vectors which 
depending on c(x) can be located anywhere in the domain (see Appendix). Since left 
singular vectors will always be to the left of their corresponding right singular vectors, left 
(right) singular vectors will always be a better approximation for B (By) than right (left). 
The norm of C~ 1 is bounded by 1 < ||£^ 1 || ? , < 2.88 from (31). So £ A is ill-conditioned 
and we expect the steady-state covariance P to have a low-rank approximation. 

We calculate P using exact dynamics and a spatial descretization with 48 grid points. 
Figure 3a shows eigenvalues Aj(P) of P along with the upper and lower bounds obtained 
from the bound matrix B and (29). Much of the variance of P is contained in a few 
modes as suggested by the ill-conditioning of £ A . The spread in the bounds is due to the 
spectrum of Po not being flat. Figure 3b shows the diagonal of P and bounds obtained 
from the bound matrix and (29). 

The relation between k and the fraction of variance explained by P^ is seen in Fig. 4; 
about 10 modes explain half of the variance. Also shown is the fraction of variance 
explained when P is projected onto the eigenvectors of B, the left singular vectors of A 
and the right singular vectors of A. The eigenvectors of B efficiently explain the variance 
of P while the left singular vectors of A do not do as well but are much better than the 
right singular vectors. 

5 Conclusions 

Ensemble and reduced-state approaches to prediction and data assimilation have shown 
low-rank covariance representations to be useful covariance models. How appropriate 
such approximation are in a given problem depends on the spectrum of the full covariance, 
which is often not available. We have examined the situation when the covariance P is a 
linear transform £ of a forcing covariance P 0 . 

The singular value decomposition of the transformation C only gives useful informa- 
tion in special cases. More useful information is obtained from operator norms of the 
transformation. Since the transformation is a positive map, mapping covariance matrices 
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to covariance matrices, there are simple expressions for its norm. Ill-conditioning is a 
necessary condition for the covariance matrix to admit a low-rank approximation. We 

i 

show that the covariance P can be ill-conditioned only when the forcing covariance P 0 
or the transformation £ are ill-conditioned. A similar result shows that the fraction of 
variance explained by the first eigenmode of the covariance can be large only when the 
same is true of the forcing covariance or when £ is ill-conditioned. 

The numerical cost of calculating these norms can be comparable to calculating the full 
covariance, though the issue of poorly known error sources is avoided. Lanczos methods 
can be used if it is possible to apply £P 0 to a vector. Theoretical insight can also be 
gained as demonstrated in the case of the discrete-time algebraic Lypunov equation where 
the usefulness of low-rank representation is directly related to stability properties of the 
dynamics. 

Appendix 

The operator A r that advances the solution of (33) in time is easily calculated by making 
a linear change of variables such that in the new variables (33) is the usual advection 
equation. The new dependent variable u(x, t ) = L (x)/j,(x, t ) where L(x) = exp (c(x)/a) 
satisfies 

u t +av x = 0, (At) 

with initial and boundary conditions 

u(x, t — 0) = L(x)fj, 0 (x), r , (x = 0,f) = 0. (A2) 

The operator A r that advances (Al) in time is A T v(x, t) = u(x — ar, t) and is related to 
A t by A t = (L(r)) _1 A T L(:r). Therefore 

A r n(x,t) = (L(x))~ l A T L(x)n(x,t) 

= L(x — ar)(L(a:)) -1 ^(a: — ar, £) (A3) 

= exp ( c(x — ar) — c(x)/a) n(x — ar, t) . 
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The singular values and vectors of A r are the eigenvalues and eigenvectors of A T A T 


A r A^(x, t) = (L(x)) 1 A(L(a:)) 2 A T (L(x)) VOM) 

= (L(x)) -1 A(L(x)) 2 (L(x + ar))~V(x + ar, t) 

(A4) 

= (L(r)) 2 (L (x ~ ar)fn(x,t) 

— exp((2c(x — ar) - 2 c(x))/a)ju(x, t ) 

except for 0 < x < ar where A r A^ is identically zero. A^A r is a multiplication operator. 
Its point spectrum contains zero and its continuous spectrum contains the range of s 2 (x ) 
for 0 < x < 1 (Halmos, 1967, Problem 66). For an eigenvalue A in the continuous 
spectrum there is no eigenfunction but there is a sequence of approximate eigenfunctions 
fj such that \\Afj — \fj\\ — > 0. Therefore, for non-constant c(x) the operator A r A^ has no 
eigenfunctions but does have approximate eigenfunctions, for example those given by the 
eigenvectors of a spatially discrete approximation. For instance, suppose the n spatial grid 
points are x* = iar/n for 0 < i < n — 1. Then, the left singular vector associated with 
the singular value s(xj) is the i-the column of the identity matrix and the right singular 
vector is the i — 1-st column of the identity matrix. 

The bound matrix B satisfies B = ABA r + I and is found by solving 


P = APA t + (L(x)) 2 
where B = L(x) _1 PL(x) _1 . P is diagonal and given by 

int(i/ar) 

p ( z ) = 51 - kar ))~ 2 > 

k = 0 

as can be verified using the relation 

in t(z/ar) 

AP(x)A t = P(x — ar) = 5^ (L(x — kar))~ 2 

k~l 


(A5) 


(A6) 


(A7) 



Therefore, 


int (x/ar) 

B(i) = (L(a:)) _2 ^ (L(x - kar )) 2 

fc =0 

int (x/ar) 

= exp —2 a~ l c(x) ^ exp 2a~ 1 c(x — kar) 

fe =0 

int(x/ar) 

= exp — 2a -1 (c(x) — c(x — kar)) 

k=0 

Similarly, the solution of B T = A* B r A + I is given by 


int(a :/ar) 

bw= exp 2a 1 {c(x) — c(x + kar)) 

k - o 


The relation 

int (x/ar) 

B(x — ar ) = exp —2a~ 1 (c(x — ar) — c(x — kar)) 

k= 1 

= B(x) exp —2 a~ l (c(x — ar) — c(x)) 

— exp — 2a~ 1 (c(:r) — c(x — ar)) 


(A8) 


(A9) 


(A10) 


shows that B(z) > B(:r — ar) if c(x) is decreasing. As a consequence the maxima of B 
and Bt occur within ar of 1 and of 0 respectively. 
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Figure Captions 


Fig. 1. c(x ) (solid line) and s(x) (dash-dot line) 

Fig. 2. (a) Mean with respect to x of the (a) forcing and (b) response plotted with respect to 
the advection time T = kra. 

Fig. 3. (a) Eigenvalues A,(P) of P (solid line) and their upper and lower bounds (dotted lines) 
obtained using (26) with Pq = A n (Po)l and Pq = Ai(Po)l. (b) Diagonal of P (solid line) and its 
bounds (dotted lines) obtained using (27) with Pq = A n (Po)l and Pq = Ai(Po)I 

Fig. 4. Fraction of the variance explained by eigenvectors of P (solid line) eigenvectors of B 
(dotted dashed line), left singular vectors of A (dashed line) and right singular vectors of A (dotted 
line). 
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(a) 



(b) 



Fig. 3. (a) Eigenvalues Aj(P) of P (solid line) and their upper and lower bounds (dotted lines) 
obtained using (26) with Pg = A n (Po)l and Pg = Ai(Po')l. (b) Diagonal of P (solid line) and its 
bounds (dotted lines) obtained using (27) with Pg = A n (P 0 )l and Pg = Ai (Pq)I 
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Fig. 4. Fraction of the variance explained by eigenvectors of P (solid line) eigenvectors of B 
(dotted dashed line), left singular vectors of A (dashed line) and right singular vectors of A (dotted 
line). 
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